function solveCaseOfStudy()
  clc
  addpath ./utils/ -end
 
  %Packet duration [slots]
  N=10;
  %Step for numeric sweep
  step=1e-3;
  %Maximum Error tolerance
  errTol=step/2;
  %Transmission rate (0-1) [packets per packet duration] 
  lambda=0.9;

  maxBOs= [ 4  3  4];               %maximum permited Backoffs Stages BOx (macMaxCSMABackoffs)
  aMinBE= [ 3  3  0];
  aMaxBE= [ 5  5  5];
  vM=     [ 4  4  4];               %vector with the number of nodes
 
  %Contention Window DO NOT CHANGE THIS VALUES
  CW=     [ 2  2  3];               

  %%Print Results
  display(['Simulation Step: ', num2str(step)]);
  display(['Error tolerance: ', num2str(errTol)]);
  display(['Class-n1 Nodes: CW=', num2str(CW(1)), ' M_bo=', num2str(CW(1)), ' macMinBE=', num2str(CW(1)), ' macMaxBE=', num2str(CW(1))]);
  display(['Class-n2 Nodes: CW=', num2str(CW(2)), ' M_bo=', num2str(CW(2)), ' macMinBE=', num2str(CW(2)), ' macMaxBE=', num2str(CW(2))]);
  display(['Class-n3 Nodes: CW=', num2str(CW(3)), ' M_bo=', num2str(CW(3)), ' macMinBE=', num2str(CW(3)), ' macMaxBE=', num2str(CW(3))]);
  
  %% Solve the Markov Model for a specific scenario, CW=2,2,3
  
  %Mean number of backoff slots waiting in BOx P^n_b
  vP{1}=vectorP(maxBOs(1), aMinBE(1), aMaxBE(1));
  vP{2}=vectorP(maxBOs(2), aMinBE(2), aMaxBE(2));
  vP{3}=vectorP(maxBOs(3), aMinBE(3), aMaxBE(3));
  
  %Find good candidates por Pi
  vPi=findPi_223(CW,maxBOs,vM,N,lambda,vP,step,errTol);
  display(['Probabilities of Channel Idle Pi=', num2str(vPi(1)), ' Pii=', num2str(vPi(2)), ' Piii=', num2str(vPi(3))]);
  
  %Solve Channel and nodes using these candidates
  [vAuxPi, vS, vPt, vPDR, pSend, vLatency, vWRx, vWTx, vWIdle, pReject]=solveModel_223(CW, maxBOs, vM, N, lambda, vP, vPi);
  display(['New probabilities of Channel Idle Pi=', num2str(vPi(1)), ' Pii=', num2str(vPi(2)), ' Piii=', num2str(vPi(3))]);
  display(['Maximum absolute error=', num2str(max(abs(vPi-vAuxPi)))]);

  %% Metrics
  display(['Aggregate Throughtput=', num2str(vS(1))]);
  display(['Throughtput               n1=', num2str(vS(2)),       ' n2=', num2str(vS(3)),       ' n3=', num2str(vS(3))]);
  display(['Prob of transmission      n1=', num2str(vPt(1)),      ' n2=', num2str(vPt(2)),      ' n3=', num2str(vPt(2))]);
  display(['Prob of sending           n1=', num2str(pSend(1)),    ' n2=', num2str(pSend(2)),    ' n3=', num2str(pSend(2))]);
  display(['Prob of Reject            n1=', num2str(pReject(1)),  ' n2=', num2str(pReject(2)),  ' n3=', num2str(pReject(2))]);
  display(['Latency                   n1=', num2str(vLatency(1)), ' n2=', num2str(vLatency(2)), ' n3=', num2str(vLatency(2))]);
  display(['Power Consumption in Transmit Mode [mW]  n1=', num2str(vWTx(1)),   ' n2=', num2str(vWTx(2)),   ' n3=', num2str(vWTx(2))]);
  display(['Power Consumption in Reception Mode [mW] n1=', num2str(vWRx(1)),   ' n2=', num2str(vWRx(2)),   ' n3=', num2str(vWRx(2))]);
  display(['Power Consumption in Idle Mode [mW]      n1=', num2str(vWIdle(1)), ' n2=', num2str(vWIdle(2)), ' n3=', num2str(vWIdle(2))]);
  display(['Total Packet Delivery Ratio=', num2str(vPDR(1))]);
  display(['Packet Delivery Ratio n1=', num2str(vPDR(2)), ' n2=', num2str(vPDR(3)), ' n3=', num2str(vPDR(3))]);


end



function mPi=findPi_223(CW,maxBOs,vM,N,lambda,vP,step,errTol)

  %Matrix with the estimated Results of pi
  mPi=[];
  tmpPi=1;
  
  error=inf;
  %Empty the Solution Vector
  mPi=[mPi;0 0 0];
  % Guessing Pi. Probability that the channel is iddle
  for pi=tmpPi:-step:step
    for piii=pi:-step:step
      %Prob that the channel is idle 'ind' consecutive slots
      %vPi(2) -> This assumption is true if and only if N1=N2=...=Nn, Eqs(5-6)
      vPi=[pi, pi*(N*pi-1+pi)/(N*pi), piii];
      %Solve the model
      [vAuxPi]=solveModel_223(CW, maxBOs, vM, N, lambda, vP, vPi);
      %Compare guessed values pi&piii with new ones computed after solving the models
      tmpErr=max(abs(vPi-vAuxPi));
      if(error>tmpErr)
        error=tmpErr;
        mPi(end,:)=vPi;
      end %if(error>
      if(error<errTol);  break; end
    end %for piii=tmpPi:
    if(error<errTol);  break; end
  end %for pi=tmpPi:
  tmpPi=mPi(end,1);
  disp(['Lambda=', num2str(lambda), ', pi=', num2str(tmpPi), ', Error=', num2str(error)]);

end

function [vAuxPi, vS, vPt, vPDR, pSend, vLatency, vWRx, vWTx, vWIdle, pReject]=solveModel_223(CW, maxBOs, vM, N, lambda, vP, vPi)
  %Maximum CW
  CWm=3;

  %Default value in case of error
  vAuxPi=ones(1,3)*NaN;

  %Probability of transmission start Eq. (1)
  p=lambda/N;

  %Cond prob. that the channel is idle, knowing that it was idle in the previous 'ind-1' backoff slots. Eqs(5-6)
  vPi_=vPi(1);
  %This assumption is true if and only if N1=N2=...=Nn
  vPi_(2)=(N*vPi(1)-1+vPi(1))/(N*vPi(1));
  vPi_(3)=vPi(3)/vPi(3-1);

  %% Solve Markov Chain for Nodes
  for i=1:CWm
    [IDLE, TX, BO, CS, tmpErr]=solveNode(p, vPi_,vP{i},maxBOs(i),CW(i));
    if(tmpErr),  return;  end
    %Summation of all possible channel states Eq(4). See Ramachandran Eq(11) for simplification details
    sumTmp=(1+TX*(N-1));
    %Probability that a Node 'ind' transmit in a generic backoff slot Eq(3)
    vPt(i)=vPi_(CW(i))*sum(CS(end,:))/sumTmp;
    if(nargout>=1)
      %Steady state probability to be in the TX state (Fraction of time spent in TX) Eq(16)
      vPtx(i)=N*TX/sumTmp;
      %Steady state probability to be in the IDLE state (Fraction of time spent in IDLE)
      vPidle(i)=IDLE/sumTmp;
      %Steady state probability to be in the BackOff, BOx, state (Fraction of time spent in BOx)
      vPbo(i)=sum(BO)/sumTmp;
      %Steady state probability to be in the Clear Channel Assessment (CCA), CSxy, state  (Fraction of time spent in CSxy)
      vPcs(i)=sum(sum(CS))/sumTmp;
      %Steady state probability to be in the first Clear Channel Assessment State (CCA), CSx1, state  (Fraction of time spent in CS1x)
      vPcs_1(i)=sum(CS(1,:))/sumTmp;

      %Prob that a node tryes to send a packet = Backoff Stages sending / Backoff stages the node has to be sending
      pSend(i)=vPtx(i)/(lambda*vPidle(i));
    end
    %Probability that a Node begin trans. given that it has sensed the channel idle 'ind' consecutive backoff slots Eq(7)
    vPt_(i)=vPt(i)/vPi(CW(i));
  end

  %% Solve Channel
  %Alpha, no one (that is allowed to) begin a transmission Eq(8)
  vA(1)=1;                                                        %no one can
  vA(2)=(1-vPt_(1))^vM(1)*(1-vPt_(2))^vM(2);                      %n1 and n2 can!!
  vA(3)=vA(2)*(1-vPt_(3))^vM(3);                                  %n1, n2 and n3 can!!

  %Beta, just one do it (mB(destination,source)) Eq(9)
  mB=zeros(3,3);

  mB(1,1)=0;                                                      %from BI1 to Sn1, cannot
  mB(1,2)=vM(1)*(1-vPt_(1))^(vM(1)-1)*vPt_(1)*(1-vPt_(2))^(vM(2));%from BI2 to Sn1, n2 can send
  mB(1,3)=mB(1,2)*(1-vPt_(3))^(vM(3));                            %from I3 to Sn1 also n3 can send

  mB(2,1)=0;                                                      %from BI1 to Sn2, cannot
  mB(2,2)=vM(2)*(1-vPt_(2))^(vM(2)-1)*vPt_(2)*(1-vPt_(1))^(vM(1));%from BI2 to Sn2, n1 can send
  mB(2,3)=mB(2,2)*(1-vPt_(3))^(vM(3));                            %from I3 to Sn2 also n3 can send

  mB(3,1)=0;                                                      %from BI1 to Sn3, cannot
  mB(3,2)=0;                                                      %from BI2 to Sn3, cannot
  mB(3,3)=vM(3)*(1-vPt_(3))^(vM(3)-1)*vPt_(3)*(1-vPt_(1))^(vM(1))*(1-vPt_(2))^(vM(2));  %from I3 to Sn3, n1 and n2 can send & from I3 to Sn2 also n3 can send
  
  [IBUSY, SUCCESS, FAIL, tmpErr]=solveChannel(vA,mB,CWm);
  if(tmpErr),  return;  end 
  
  %% Compute vAuxPi from channel
  %Sum of dwell time in all channel states Eq(13)
  sumChSt =sum(IBUSY)+N*(sum(SUCCESS)+FAIL);

  %Cond prob that the channel is in the Busy|Idle_i state Eq(12)
  vPi_b=IBUSY(1:CWm-1)/sumChSt;

  %Prob that the channel is idle CW slots Eq(12)
  vAuxPi(CWm)=IBUSY(end)/sumChSt;

  %Cond prob. that the channel is idle, knowing that it was idle in the previous 'ind-1' backoff slots Eq(15)
  for i=1:CWm-1
    %Prob that the channel is idle 'ind' consecutive slots
    vAuxPi(i)=vAuxPi(CWm)+sum(vPi_b(i:CWm-1));
  end
  
  %% Just vAuxPi??
  if(nargout==1)
    return;
  end

  %% Results
  display(['     Transition Probabilities Alpha= ', num2str(vA)]);
  display(['     Transition Probabilities to a Success state Beta-1= ', num2str(mB(1,:))]);
  display(['                                                 Beta-2= ', num2str(mB(2,:))]);
  display(['                                                 Beta-3= ', num2str(mB(3,:))]);
  
  %Throughput [Stot, S1, S2, ..., Sn] per class. Effective speed (KBps) = Sx/Mx* 250e3Kbps/(8b/B) Eq(16)
  vS=N*SUCCESS/sumChSt;
  vS=[sum(vS), vS];

  % Packet Delivery Ratio Type 2: #successfully received packets/#sent packets = sum(S)/sum(Mi*Ptxi) Eq(26)
  vPDR=[vS(1)/sum(vPtx.*vM), vS(2:end)./(vPtx.*vM)];

  %Latency in backoff slots Eq(27)
  vLatency=N*vM./vS(2:end).*(1-vPidle);

  %Power consumption Chipcon cc2420 802.15.4 compliant RF transceiver
  Nbeacon=2;              % Beacon length (in backoff slots)
  Nir=0.6;                % This is the duration of the switching period idle->receive (in backoff slots)
  BI=3072;                % Beacon Interval (in backoff slots)
  Widle=712e-3;           % Power consumed in idle state (mW)
  Wrx=35.28;              % Power consumed in RX state (mW)
  Wtx=31.32;              % Power consumed in TX state (mW)

  Pbeacon=Nbeacon/BI;     %Probability to receive a beacon -> Beacon Length/Beacon Interval length
  Pir=Nir*(1/BI+vPcs_1);  %Prob to be activating the radio Time to switch idle-receive/BI due to be in CSx1 or reciving Beacon

  %Power cosumed during reception Eq(20)
  vWRx=(vPcs+Pbeacon+Pir)*Wrx;
  %Power cosumed during transmission Eq(19)
  vWTx=vPtx*Wtx;
  %Power cosumed during idle period Eq(18)
  vWIdle=(vPidle-Pbeacon-Pir+vPbo)*Widle;

  %Probability of throw away a new packet Eq(23)
  pReject=(1-vPidle);  
end

